#include <stdio.h>
#include "splines.c"
#include "mpi.h"

#include "randomlib.h"
#include "randomlib.c"
#include "fminbr.c"
#include "time.h"

/* double fminbr(double a,double b, double (*f)(double x),double tol); */
double f(double x);
void dat_generate(struct L2_1D_DATA*, int, int, double);
void gnomesort(int n, double * ar);
void get_random_knots(struct SP_1D * sp, struct L2_1D_DATA * data );
double get_L2_error(struct SP_1D * sp, struct L2_1D_DATA * data );
double opt_knot(void);
struct SP_1D * par_fit(int seed, int order, int n, double * glob_sse);
struct SP_1D * get_one_less_knot(int i, struct SP_1D * sp);
struct SP_1D * step_down_fit( int order, int n_large, int n_final, double * glob_sse);
struct SP_1D * step_down_iter(int num_tries, int order, int n_large, int n_final);
double fun(double x, int i);


/* Global data for each process*/
int iknot;
struct L2_1D_DATA *  dat;
struct SP_1D * sp, * save_sp[1000];
int  num_rand_splines;


int main(int argc, char* argv[])
{


  int my_rank;     /* rank of process */
  int p;           /* number of processes */
  
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &p);

  /* Declare the global array for storing all the output data */
  double *global_spec_loss;
  double *global_maic_loss;
  int *global_maic_order;
  int *global_maic_nknots;

  /* 
    ------------------------------------------
     Starting the instructions on each process 
    ------------------------------------------- 
  */

  FILE *fp, *fp1;

  /* n_final is the final dimension of the spline, when n_final equals order then there is no interior knots */
  int i, k, idx, imax, order, order_min, order_max, order_spec, nknots_spec, n, nobs, inform, n_large, n_final, off_set, order_vec[500], rank[500], nknots_vec[500];
  int mod_aic_min_order[100], mod_aic_min_nknots[100];
  int seed1, seed2, count, sum, i_maic, i_maic_spec, i_true_ssq, num_step_dn, num_fn, num_dat;


  struct SP_1D  * sp1;
  struct L2_1D_DATA *  dat_actual;

  double h, x, y, sse, cv, cor, glob_sse=999999, glob_sse_loc=999999., sigma, act_err, aic, bic;
  double mod_aic, mod_aic_vec[500], act_err_vec[500], sse_vec[500],  dim_vec[500], sort_v[500];
  double min_maic, min_true_ssq, min_maic_spec, min_loss_spec;
  double mod_aic_min[100], act_err_min[100], mod_aic_spec[100], act_err_spec[100];
  float read_sigma, ratio;
  clock_t t1,t2;


  /* Start the clock for the root processor. */
  if(my_rank == 0){
    t1 = clock();
  }
  /* Open the files. */
  fp = fopen("../spa6/step_output_spa6_800.txt", "w");
  fp1 = fopen("../spa6/step_input_spa6_800.txt", "r");

  /* fprintf(fp, " My rank is %d", my_rank); */

  /* Read the input parameters. */
  fscanf (fp1, " %d %d %f",  &num_fn, &nobs, &read_sigma);
  fscanf (fp1, "%d %d %d", &count, &num_step_dn, &n_large);
  fscanf (fp1, "%d %d", &order_spec, &nknots_spec);
  fscanf (fp1, "%d %d %d", &order_min, &order_max, &imax);

  /* Initialize some parameters. */
  sigma = read_sigma;
  h = 1./(nobs -1);
  order = 0;
  ratio = 1./CLOCKS_PER_SEC;

  /* Declare the local loss arrays. */
  double local_spec_loss[imax];
  double local_maic_loss[imax];
  int local_maic_order[imax];
  int local_maic_nknots[imax];
  
  /* Initialize the data_actural. */
  dat_actual = data_1d_initialize(  nobs, &inform);
  dat_generate(dat_actual, 0, num_fn, 0);


  /* -------------------------------------- 
     Start the major loop for each process. 
     --------------------------------------*/

  /* The level-1 loop iterates the seeds*/
  for(idx=0; idx<imax; idx++){ 

    /* Initialize the data struct. */
    dat = data_1d_initialize(  nobs, &inform) ;
    num_dat = idx + my_rank*imax;

    /* Generate the data struct */
    dat_generate(dat, num_dat, num_fn, sigma);

    /* k is increased by 1 for every inner loop iteration */
    k= -1;    

    /* The level-2 loop iterates orders*/
    for(order = order_min; order <order_max+1; order++){
      /* The level-3 loop iterates dimensions */
      for(n_final = order; n_final < count+order; n_final++){
	k++;
	/* 
	   step_down_iter is the core function that go through 
	   the step_down phase and opt_knot phase.

	  Starting dimension: n_large;
	     Starting nknots: n_large-order;
	   Ending dimenstion: n_final;
	       Ending nknots: n_final-order;
	*/
	sp1 = step_down_iter(num_step_dn,  order, n_large, n_final);

	sse = get_L2_error( sp1, dat);

	save_sp[k] = spline_1d_copy(sp1);                // Store the spline

	act_err = get_L2_error( sp1, dat_actual);
	aic = dat->nobs*log(pow(sse,2)/dat->nobs) + 2*(2*sp1->n - sp1->order);
	bic = dat->nobs*log(pow(sse,2)/dat->nobs) + log(dat->nobs)*(2*sp1->n - sp1->order);
	mod_aic = aic + (2*sp1->n - sp1->order);

	mod_aic_vec[k] = mod_aic;                        // Store the maic
	act_err_vec[k] = act_err*act_err/dat->nobs;      // Store the loss
	sse_vec[k] = sse*sse/dat->nobs;                  // Store the error
	order_vec[k] = sp1->order;                       // Store the order
	nknots_vec[k] = n_final - sp1->order;            // Store the nknots

	free_SP_1D(sp1);
      }
    }

    i_maic = 0;
    i_maic_spec = 0;
    i_true_ssq = 0;
    min_maic = mod_aic_vec[0];
    min_true_ssq = act_err_vec[0];

    /* Search for the minimum MAIC and minimu loss over all orders. */
    for(i=1; i< count*(order_max-order_min+1); i++){
      if(min_maic > mod_aic_vec[i]){
	i_maic = i;                             // Store the min maic index
	min_maic =  mod_aic_vec[i];             // Store the min maic value
      }
      if(min_true_ssq > act_err_vec[i]){
	i_true_ssq = i;                         // Store the min loss index
	min_true_ssq =  act_err_vec[i];         // Store the min loss value
      }
    }
    
    /* Complute the loss for specific order and nknots. */
    min_loss_spec = act_err_vec[count*(order_spec-2)+nknots_spec];

    /* Search for the minimum MAIC over specific order. */
    min_maic_spec = mod_aic_vec[0];
    for(i=count*(order_spec-order_min); i<count*(order_spec-order_min+1); i++){
      if(min_maic_spec > mod_aic_vec[i]){
	i_maic_spec = i;                        // Store the specific maic index
	min_maic_spec =  mod_aic_vec[i];        // Store the specific maic value
      }
    }

    /* Search for the rank k */
    for(i=0; i< count*(order_max-order_min+1); i++)
      sort_v[i] = act_err_vec[i];    
    gnomesort((order_max-order_min+1)*count, sort_v);
    k = 0;
    for(i=0; i<(order_max-order_min+1)*count; i++){
      if (sort_v[i] == act_err_vec[i_maic]) break;
      k = k++;
    }


    /* Begin output  */
    /* fprintf(fp, "\n \n M_aic spline has rank = %d \n \n", k); */

    /* Echo the input parameters  */
    /*
      fprintf(fp, "\n \n The rn parameters are: \n");
      fprintf(fp, "seed1 = %d, seed2 = %d, sigma =  %f \n \n", seed1, seed2, sigma);
      fprintf(fp, "Function number = %d  and the value at 0 is  %f \n", num_fn, fun(0., num_fn));
      fprintf(fp, "NOBS = %d \n \n", dat->nobs);
      fprintf(fp, "order = %d \n \n", sp1->order);
      fprintf(fp, "num_step_dn = %d ", num_step_dn);
      fprintf(fp, "n_large = %d \n \n \n", n_large);

      if(i_true_ssq == i_maic)
      fprintf(fp, "M_aic correcly predicts min IMS \n \n");
      else
      fprintf(fp, "No agreement between best fit and aic \n \n");

      fprintf(fp, "act_err/nobs  mod_aic      SSE/nobs  \n");

      i = 0;
      for(order = 2; order <7; order++){
      for(k= 0; k< count; k++){
      fprintf(fp, "%e     %f   %f  %d\n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);
      i++;}
      fprintf(fp, "\n");
	
      }
      sum = 0;
      i = 0;
      for(order = 2; order <7; order++){
      for(k=0; k<count-2; k++){
      if(sse_vec[i+1] > sse_vec[i])
      {sum++; break;}
      i++;}
      if (sum > 0 ) break;
      i = count*(order-1);
      }
	
      if(sum >0) {
      fprintf(fp, "                           Caution errors NOT ordered  \n");
      fprintf(fp, "                           order = %d   position = %d  %d \n", order, k, i);
      }
      else
      fprintf(fp, "                           Errors ordered as expected \n");
      fprintf(fp, "\n \n");
      fprintf(fp, "act_err/nobs  mod_aic      SSE/nobs    Spline order\n \n");
      k = sp->order;
      i = i_true_ssq;
      fprintf(fp, "Actual Best Error \n");
      fprintf(fp, "%f     %f      %f         %d \n \n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);
      i = i_maic;
      fprintf(fp, "Minimum Modified AIC \n");
      fprintf(fp, "%f     %f      %f         %d \n \n \n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);
    */

    /* Record the best and special fit into vectors */
    act_err_spec[idx] = min_loss_spec;                   // Store the specific loss
    mod_aic_min[idx] = min_maic;                        // Store the min maic 
    act_err_min[idx] = act_err_vec[i_maic];             // Store the loss of the min maic 
    mod_aic_min_order[idx] = order_vec[i_maic];                 // Store the order of min maic spline
    mod_aic_min_nknots[idx] = nknots_vec[i_maic];               // Store the nknots of min maic spline

    /* mod_aic_spec[idx] = mod_aic_vec[i_maic_spec]; */
    /* act_err_spec[idx] = act_err_vec[i_maic_spec]; */

  } /* End of major loop */

  /* Record the end time and print the elapsed time in seconds. */
  
  /* 
     -----------------------------------
     End of instructions on each process.
     ----------------------------------- 
  */



  /* Fill out the local loss arrays. */
  for(i=0; i<imax ; i++){
    local_spec_loss[i] = act_err_spec[i];
    local_maic_loss[i] = act_err_min[i];
    local_maic_order[i] = mod_aic_min_order[i];
    local_maic_nknots[i] = mod_aic_min_nknots[i];
  }

  /* Allocate space for the global arrays. */
  if(my_rank == 0){
    global_spec_loss = (double*)malloc(imax*p*sizeof(double));
    global_maic_loss = (double*)malloc(imax*p*sizeof(double));
    global_maic_order = (int*)malloc(imax*p*sizeof(int));
    global_maic_nknots = (int*)malloc(imax*p*sizeof(int));
  }

  /* Gather the local loss arrays to fill the global loss arrays. */
  MPI_Gather(local_spec_loss, imax, MPI_DOUBLE, global_spec_loss, imax, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Gather(local_maic_loss, imax, MPI_DOUBLE, global_maic_loss, imax, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Gather(local_maic_order, imax, MPI_INT, global_maic_order, imax, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Gather(local_maic_nknots, imax, MPI_INT, global_maic_nknots, imax, MPI_INT, 0, MPI_COMM_WORLD);


  /* Print the global loss arrays. */
  if(my_rank == 0){
    t2 = clock();
    fprintf(fp, "\n Running Time = %f \n", ratio*(long)t1+ratio*(long)t2);
    fprintf(fp, "The specific loss        The loss by MAIC\n\n");
    for(i=0; i<imax*p; i++)
      fprintf(fp, "%f %f %d %d\n", global_spec_loss[i], global_maic_loss[i], global_maic_order[i], global_maic_nknots[i]); 
    free(global_spec_loss);
    free(global_maic_loss);
    free(global_maic_order);
    free(global_maic_nknots);
  }

  MPI_Finalize();

  return 0;
}


void dat_generate(struct L2_1D_DATA* dat, int num_data, int num_fn, double sigma){

  int i;
  double x;
  double h = 1./(dat->nobs - 1);

  int seed1 = 100 + 10*num_data;
  int seed2 = 1000 + 5*num_data;

  RandomInitialise(seed1, seed2);

  for(i=0; i<dat->nobs; i++){
    dat->weight[i] = 1;
    dat->xdata[i] = i*h;
    x = dat->xdata[i];
    dat->ydata[i] = fun(x, num_fn) + sigma*RandomGaussian(0., 1.);
  }
}


void get_random_knots(struct SP_1D * sp, struct L2_1D_DATA * data ){
  int i, num_int_knots;
  double *  temp, low, high;

  num_int_knots = sp->n - sp->order;
  temp = (sp->knot + sp->order);

  /*  Note this assumes that the xdata field is ordered */
  low = data->xdata[0];
  high = data->xdata[data->nobs -1];

  for(i=0; i<num_int_knots; i++)
    temp[i] = low +RandomUniform()*(high - low);
	
  gnomesort(num_int_knots, temp);

  /*
    for(i=0; i<num_int_knots; i++)
    sp->knot[sp->order + i] = temp[i];

    free(temp);
  */
}


void gnomesort(int n, double * ar){
  int i = 0;
  while(i<n){
    if (i==0 || ar[i-1] <= ar[i]) i++;
    else {double tmp = ar[i]; ar[i] = ar[i-1]; ar[--i] = tmp;}
  }

}

double get_L2_error(struct SP_1D * sp, struct L2_1D_DATA * data ){
  int i;
  double sse;	
  sse = 0.;
  for(i=0; i<data->nobs; i++)
    sse += (data->weight[i])*pow(data->ydata[i] - sp_1d_value(sp, data->xdata[i], 0), 2);

  return sqrt(sse);
}

double opt_knot(void){
  int   j;
  double sse_old, sse, a, b, tol; 

  /* tol = .0001; */
  tol = .0000001;
	
  sse_old = get_L2_error( sp, dat );
  for(j=0; j<99; j++){
    for(iknot = sp->order; iknot < sp->n  ; iknot++){
      a = sp->knot[iknot-1], b = sp->knot[iknot+1];
      sp->knot[iknot] = fminbr(a, b, f, tol);
	
    }

    for(iknot = sp->n -1; iknot >= sp->order  ; iknot--){
      a = sp->knot[iknot-1], b = sp->knot[iknot+1];
      sp->knot[iknot] = fminbr(a, b, f, tol);
		
    }
	

    for(iknot = sp->n -1; iknot >= sp->order  ; iknot--){
      a = sp->knot[iknot-1], b = sp->knot[iknot+1];
      sp->knot[iknot] = fminbr(a, b, f, tol);
		
    }

    for(iknot = sp->order; iknot < sp->n  ; iknot++){
      a = sp->knot[iknot-1], b = sp->knot[iknot+1];
      sp->knot[iknot] = fminbr(a, b, f, tol);
	
    }
    sse = get_L2_error( sp, dat );
    /*	printf( " sse = %f   sse_old = %f \n", sse, sse_old);  */
    if( fabs( (sse-sse_old)/(sse_old + tol) ) < tol )
      break;
    sse_old = sse;
  }

  sse = get_L2_error( sp, dat );
  return sse;
}

double f(double t){
  double sse;
  sp->knot[iknot] = t;
  get_L2_1D_spline(dat, sp);
  sse = get_L2_error( sp, dat );
  return sse;
}

	

struct SP_1D * par_fit(int seed, int order, int n, double * glob_sse){
  int i,j, k, inform;

  struct SP_1D  * sp1, *sp2;

  double   sse, glob_sse_loc=999999.;

  /* Data and sp are global structures */

  *glob_sse=999999;


  k = seed;

  /*
    printf(" \n \n The seeds for this job are %d and 15000\n", k);
    printf("The number of random splines is %d per iteration \n", num_rand_splines);
    printf("Order = %d and dimension = %d \n\n ", order, n);

  */

  sp = sp_1d_initialize( order, n,  &inform);

  sp_1d_set_knots(sp, 0., 1.);

  get_L2_1D_spline( dat, sp );
  sp1 = spline_1d_copy(sp);
  sp2 = spline_1d_copy(sp);


  RandomInitialise(k,15000);

  /* Begin major loop  */

  for(k=0; k<9; k++){

    glob_sse_loc=999999.;

    for(j=0; j< num_rand_splines; j++){
      get_random_knots(sp, dat);



      get_L2_1D_spline( dat, sp );

      sse = get_L2_error( sp, dat);





      if(sse < glob_sse_loc){
	free_SP_1D(sp2);
	glob_sse_loc = sse;
	sp2 = spline_1d_copy(sp);}
		
      /*	  printf("The error is %f \n", sse);  */

    }


    /*
      printf("The iteration is %d \n", k);

      printf("The global best error before opt is %f \n", glob_sse_loc);
      for(i=0; i< sp2->n + sp2->order; i++)
      printf("%f ", sp2->knot[i]);

      printf("\n ");

    */

    free_SP_1D(sp);
    sp = spline_1d_copy(sp2);

    sse = opt_knot(); 

    if(sse < *glob_sse){
      free_SP_1D(sp1);
      *glob_sse = sse;
      sp1 = spline_1d_copy(sp);
	
      /*
	printf("The global best error after opt is %f \n", *glob_sse);
	for(i=0; i< sp1->n + sp1->order; i++)
	printf("%f ", sp1->knot[i]);


      */


    }

    printf("\n \n");

  }

  return   sp1;
}




struct SP_1D * step_down_fit( int order, int n_large, int n_final, double * glob_sse){
  int i,j, k, isave, inform;

  struct SP_1D  * sp1, *sp2;

  double   sse, glob_sse_loc=999999.;

  /* Data and sp are global structures */

  /*
    printf(" \n \n The seeds for this job are %d and 15000\n", k);
    printf("The number of random splines is %d per iteration \n", num_rand_splines);
    printf("Order = %d and dimension = %d \n\n ", order, n);

  */

  sp = sp_1d_initialize( order, n_large,  &inform);

  sp_1d_set_knots(sp, 0., 1.);

  get_L2_1D_spline( dat, sp );
  sp1 = spline_1d_copy(sp);
  sp2 = spline_1d_copy(sp);
  get_random_knots(sp,  dat );

  /* Begin major loop  */

  k = n_large;

  while( k > n_final){

    glob_sse_loc=999999.;
	
    /*	if(5*(k/5) == k)
	sse = opt_knot(); 
    */
    for(i = order; i<k; i++){

      free_SP_1D(sp2);
      sp2 = get_one_less_knot(i, sp);	
      get_L2_1D_spline( dat, sp2 );

      sse = get_L2_error( sp2, dat);


      if(sse < glob_sse_loc){
	isave = i;
	free_SP_1D(sp1);
	glob_sse_loc = sse;
	sp1 = spline_1d_copy(sp2);}
		
      /*	  printf("The error is %f \n", sse);  */

    }
    //printf("For k = %d  The error is %f \n", k, glob_sse_loc);
    k--;
    free_SP_1D(sp);
    sp = spline_1d_copy(sp1);

  } //end while
  /*
    printf("The iteration is %d \n", k);

    printf("The global best error before opt is %f \n", glob_sse_loc);
    for(i=0; i< sp2->n + sp2->order; i++)
    printf("%f ", sp2->knot[i]);

    printf("\n ");

  */

  free_SP_1D(sp);
  sp = spline_1d_copy(sp1);

  sse = opt_knot(); 

  if(sse < *glob_sse){
    free_SP_1D(sp1);
    *glob_sse = sse;
    sp1 = spline_1d_copy(sp);	}
	
  /*
    printf("The global best error after opt is %f \n", *glob_sse);
    for(i=0; i< sp1->n + sp1->order; i++)
    printf("%f ", sp1->knot[i]);


  */




  //printf("\n \n");
  return   sp1;

}





struct SP_1D * get_one_less_knot(int i, struct SP_1D * sp){
  struct SP_1D * sp1;
  int j, inform;
	
  sp1 = sp_1d_initialize( sp->order, sp->n -1,  &inform);

  for(j=0; j< sp->n + sp->order; j++){
    if(j < i)
      sp1->knot[j] = sp->knot[j];
		
    else if(j>i)
      sp1->knot[j-1] = sp->knot[j];
		
  }

  return sp1;	
}




struct SP_1D * step_down_iter(int num_tries, int order, int n_large, int n_final){
  int i, j;
  double glob_sse = 999999, cur_sse;
  struct SP_1D * sp_final, * sp1;

  cur_sse = glob_sse;
  for(i=0; i< num_tries; i++){
    sp1 = step_down_fit(order, n_large, n_final, &glob_sse);
    if(glob_sse < cur_sse){
      sp_final = spline_1d_copy(sp1); 
      cur_sse = glob_sse;} 

  }
  free_SP_1D(sp1);
  return sp_final;
}

double fun(double x, int i){
  const double pi = acos(-1);  
  double y= 0.;
  switch(i){
  case 1: y = 10*(sin(x*x*x*12  ));  break;
  case 2: y = 1/(.1 + sin(x*x*8 )*sin(x*x*8 )); break;
  case 3: y = fabs(cos(x*x*6)); break;
  case 4: y = 3*exp( -50*(x-.1)*(x-.1)) + 6*exp( -100*(x-.7)*(x-.7));break;
  case 5: y = 10*(pow(fabs(x-.2), 3) - pow(fabs(x-.3), 3) + 4*pow(fabs(x-.5), 3) - pow(fabs(x-.8), 3)); break;
  case 6: y = (fabs(x-.2) + fabs(x-.3) + fabs(x-.5) + fabs(x-.8)); break;
  case 7: y = sin((1./(x*x+.1)) + 1./(x-1.2)); break;
  case 8: y = 8*exp( -100*(x-.1)*(x-.1)) + 
      6*exp( -200*(x-.6)*(x-.6))+ 
      10*exp( -150*(x-.9)*(x-.9));break;
  case 9: y = 2*sin(8*exp( -100*(x-.1)*(x-.1))) + 
      6*exp( -200*(x-.6)*(x-.6))+ 
      10*exp( -150*(x-.9)*(x-.9));break;
    /* SpaHet3 */
  case 10: y = sqrt(x*(1-x))*sin((2*pi*(1+pow(2,-3/5)))/(x+pow(2,-3/5)));break;
    /* SpaHet6 */
  case 11: y = sqrt(x*(1-x))*sin((2*pi*(1+pow(2,-3)))/(x+pow(2,-3)));break;
  }
  return y;
}
